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Abstract 

Exchange symmetry in acceleration partitions the configuration space of an N particle, one- 
dimensional, gravitational system into AM equivalent cells. We take advantage of the resulting 
small angular extent of each cell to construct a related integrable version of the system that takes 
the form of a central force problem in iV — 1 dimensions. The properties of the latter, including 
the construction of trajectories, as well as several continuum limits, are developed. Dynamical 
simulation is employed to compare the two models. For a class of initial conditions, excellent 
agreement is observed. 
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One dimensional systems have played a central role in the development of physical models. 
Because of their relative simplicity, they form the starting point for the study of a wide 
variety of interactions. One-dimensional models of gases have been studied extensively in 
statistical mechanics and one- dimensional lattice models have been used in solid state physics 
to study metal alloys, magnetic spin systems, glasses, phonon propagation, electronic bands 
and a host of other systems. One dimensional gravitational systems have been used to 
investigate relaxation to equilibrium, phase transitions, the gravothermal catastrophe, and 
cosmology Q]. An important example is the system consisting of a one- dimensional chain 
of coupled, nonlinear, oscillators. Known as the Fermi-Pasta-Ulam (FPU) problem^], it 
has had an enormous impact on the development of dynamics in the last fifty years. Its 
failure to come to thermal equilibrium in simulations on the earliest electronic computers 
stimulated intense research in the theory of conservative dynamical systems and resulted 
in major breakthroughs such as the KAM theorem and developments in soliton theory Q]. 
The analogous system in astronomy consists of iV parallel planar mass sheets of infinite 
extent that are restricted to move perpendicular to their surface. Like the FPU problem, 
this is a nonlinear, one-dimensional system. The only forces on a given sheet (or particle) 
are gravitational, and the system obeys the Poisson equation in one dimension. Since each 
particle experiences a constant acceleration until it encounters a neighboring particle, the 
motion can be expressed analytically and it is straightforward to simulate the system on the 
computer jsj. 

The one dimensional self-gravitating system (OGS) was the first gravitational system to 
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equilibrium. Early studies showed that the system virialized within 50 or so characteristic 
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system crossing times, r c , as a consequence of "violent relaxation" [5J. Following this 

period the system appears to approximate a stationary Vlasov state which evolves extremely 
slowly, perhaps on the order of millions of r c . Evidence of the approach to thermal equilib- 
rium has only been recently established for a two component system . For an extensive 
review of the history of the relevant investigations see L3| . In this paper we offer a possible 
explanation for the incredibly slow exploration of phase space these systems exhibit. We 
demonstrate analytically and with numerical simulation that an exactly integrable Hamil- 
tonian system (EIS) exists which closely represents the OGS as its population increases. 
Further we derive a continuum (Vlasov) limit for the EIS and show that it is characterized 
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by two integral invariants, as well as an exactly periodic radius. We also show that the 
continuum model can be extended to the pair space and provide information concerning 
correlations in position and velocity. 

Consider an OGS consisting of N parallel mass sheets (particles) of equal mass surface 
density m and let the real x axis be perpendicular to their surface. Label the sheets in terms 
of their increasing ordering from left to right so their positions are Xj with Xj + i > Xj(j = 
1,...,N — 1) . Then, from Gauss' law, the acceleration Aj experienced by particle j is a 
constant proportional to the net difference in the number of particles to its right and left: 

A j = 2irmG(N+l-2j). (1) 

As the system evolves, each particle undergoes uniform acceleration until a pair of par- 
ticles crosses, at which time their accelerations are exchanged. Since the system is isolated 
from external forces, the total momentum and energy are conserved, while the center of 
mass moves with uniform velocity. Without loss of generality we may assume we are in a 
frame where the total momentum is zero and the center of mass is at rest. To establish a 
connection with the Vlasov limit, it is more convenient to represent the system in terms of 
positions and velocities (xj,Vj) of the particles. Thus 

N N N 

5>* = o> E^' = > £ = yE^+ WG Ei^-^i ( 2 ) 
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where E is the total energy. It is also convenient to introduce a system of unit vectors , 
ej , in the N dimensional configuration space. Then we can represent the system in terms 
of N dimensional position, velocity (in the tangent space), and acceleration vectors, 

N N N 

x = Y x i e 3 ' v = Y v i e i ' A = Y A i e r ( 3 ) 
ill 

In this context, the system is represented by a single point moving with constant accel- 
eration until a crossing occurs, whereby two of the acceleration components are exchanged; 
i.e. the actual components of A depend on the ordering of the particle positions. Letting 
e c =-L Ylj e j> the constraints on the total momentum and the center of mass can be expressed 
by v ■ e c = x ■ e c = 0. As a result of the dynamical constraints (including energy conserva- 
tion), the trajectories lie on a 2N — 3 hyper-surface of the 2N dimensional phase space. We 
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can also introduce a generalized angular momentum tensor as a two form, or dyadic in the 
older formalism |8j, and its square: 



N N 

L = x A mv = m XiVj&i A Bj = —= (xiVj — XjVi)eiej (4) 

i,j=l V i,j=l 

2 N 

L 2 = — ( X i V i ~ X 3 V i) 2 - ( 5 ) 

It is straightforward to show that, for iV = 3, these definitions reduce to the familiar 
equations of three dimensional orbital mechanics. For completeness, and to compare with 
our further development below, we mention that the system has a well defined continuum 
(Vlasov) limit. By projecting the single point representing the system in its 2N dimensional 
phase space onto a single (x, v) plane, the system is represented by a cloud of N points. If 
we let N — > oo while constraining the total mass, M = Nm , and energy E, the system 
is described by a fluid in the f-i(x,v) space with normalized distribution function f(x,v;t) 
satisfying the Vlasov (or Vlasov-Poisson) equation 

df df . ,df 

^ + ^ + A[ % = °' (6) 
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A[f] =2nMG J dv'J dx'f(x',v',t)[Q(x' - x) - Q(x - x')] (7) 

where 6(x)is the usual step function jjllll- This system has been studied in t 
extensively and is known to have an infinite number of stationary solutions 
particular, the maximum entropy solution is / ~ exp (-/3t» 2 /2) cosh- 2 (27rM 2 Gx/2£). 

To understand our approach, we first consider a simpler system of a single particle moving 
in the plane in which we introduce the usual polar coordinates, (r, We imagine that 
the plane is segmented into P pie shaped "wedges" which share a common vertex at the 
origin and subtend 2tt/P radians. In the k th wedge the particle experiences a constant 
acceleration g^, with |gfc| = g,(k = 1, P) , directed along the wedge bisector, thus pointing 
approximately towards the origin (common vertex, see Fig. Then within each segment 
the particle executes a parabolic trajectory until it reaches a wedge boundary, which it 
crosses and begins a new parabolic segment with a rotated acceleration. If, as in Fig. 
P = 6 it has been shown that this system is isomorphic to the system of three parallel mass 
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sheets Q|. Note that the constraint on the center of mass constrains the motion to the 
plane. We now imagine what happens as P becomes large. In each wedge the acceleration 
is directed more nearly in the radial direction, i.e. towards the origin. Thus the system is 
very close to another one where the acceleration g = —gr/r varies smoothly and is defined 
everywhere except at the origin. This latter system is exactly integrable: both energy and 
angular momentum are exactly conserved. For large, but finite, P, most orbits in the former 
system will lie close to their integrable counterparts for long times. However, for any P < oo, 
there will be time intervals where % changes sign for some of the orbits in the segmented 
system and the orbits will start to diverge from their integrable cousins. 

Now let's shift our attention back to the OGS. Note that there are N\ possible orderings of 
the particle labels. Associated with each ordering is a pyramidal segment in the configuration 
space in which the acceleration, A, is a unique constant. To explore how the geometry 
changes with increasing N, let's compute cos(0) := A • A* /A 2 where A* has the permuted 
particle ordering to A obtained by interchanging the positions of a pair of adjacent particles, 
and A 2 = ^2 A 2 . With a little work we obtain, 

2{2 *f A f = sm 2 (0/2) , A 2 = \(2,MG) 2 N{1 - AT 2 ). (8) 

Therefore, for N ^> 1, = 2\^6/N 3 ^ 2 becomes vanishingly small, and we are faced with 
a similar situation to the planar example discussed above, i.e. there are N\ cells within 
which the acceleration is approximately in the radial direction in the N — 1 dimensional 
space obtained by projecting out the center of mass direction e c . As a consequence we 
can identify the OGS with a different system consisting of a single particle moving in N 
dimensions with a radial acceleration of constant magnitude A. As we shall quickly see, the 
latter is an exactly integrable system. Let us signify the coordinates and velocities of the 
smooth system by (yj,Uj) so its position and velocity are y =J2Vj e j^ u = J2 u j e j- Then 
the new equations of motion are 



dyj duj Ayj 



dt ' 1 ' dt ' R 



N 



This is a central force problem in N dimensions with a potential energy <3> = mAR. Direct 
substitution in the above quickly shows that the energy, H := |k 2 + $, the generalized 
angular momentum two form L =my A v, and L 2 are exact dynamical invariants of the 
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FIG. 1: Partition of the plane into six segments. Arrows indicate direction of acceleration. Note 
that within each wedge shaped segment the acceleration is pointing inward and exactly parallel to 
the wedge bisector. 

motion. Moreover, if at an initial time, e c -y = e c u = 0, then by taking successive derivatives 
of y and u with respect to t it can be seen that these quantities don't change their value, 
and the trajectories remain constrained to an iV — 1 dimensional manifold. It is remarkable 
that, in spite of the arbitrary number of dimensions, since % ~ y , the motion remains in 
the two-plane denned by the initial vectors yo and u ! Therefore the motion is simple to 
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describe: If we impose polar coordinates (R, 6) in the plane of the motion we obtain 



d 2 R L 2 d9 _ L 



dt 2 m 2 R 3 ' dt mR 2 

Alternatively, a first order equation for R(t), 



dR 
~dt 




L 2 

H - mAR - 



(11) 



2mR 2 

can be obtained directly from the energy equation and shows that turning points of the 
orbit occur where the argument of the radical vanishes. Consequently R(t) is a periodic 
function of the time with period, say, t r . However, we must be careful to note that 6{t) 
does not advance by 2n radians in the time Tr , so the orbits are not closed in space. By 
inverting the equation (solving for t = t(R)) it can be shown that the solution is a linear 
combination of elliptic integrals. 

To obtain the complete set of coordinates and velocities, we use the initial position and 
velocity, y = y(t = 0), u = u(t = 0), to define orthogonal unit vectors in the plane, bi , 
b 2 : 

u y° u _ (up - (up ■ e y )e y ) 
R (tig - (u • e^) 2 ) 

Then we can impose Cartesian coordinates (z±, z 2 ) such that 

y(t) = zibi + z 2 h 2 , u(t) = ^bi + ^b 2 (13) 

where {ziz 2 ) satisfy 

d2zi - a Zl d2z2 - a Z2 ( U ) 

d ^ ~ ' dt2 " 

The complete set of coordinates and velocities are obtained from yj(t) = ej ■ y(t) , 

Uj(t) = ej ■ u(t). Note that the Cartesian pair (zi : z 2 ) can be obtained either by solving this 

pair of coupled second order equations, or by simply noting that z\ = Rcos9, z 2 = Rsin9. 

Several continuum limits for the exactly integrable gravitational analogue (EIS) can be 

constructed by focusing on R 2 =^2 y 2 . As for the OGS, project the phase space onto a single 

two dimensional plane (y, u) to obtain a cloud of N points, and take the limit iV — ► oo 

while constraining the total mass, M = Nm , and energy H. Let fi(y, u; t) be the normalized 
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density pdf in the (y, u) plane. Then in the continuum limit we can immediately identify 
R 2 = YliV'j w ith N J dy J duy 2 fi := N (y 2 ) which, in general, depends on the time. We will 
see that in the large N limit the factors of N will scale out of the equations. Recalling that 
A ~ \JlSf , the equations of motion of a point in the (y, u) plane in the continuum limit are 

dy du Ay (27iMG)y 

~r — u i —r — = , 15) 

dt dt ^/ww) VW) 

which results in the following Vlasov equation for fi(y, u, t): 

dj, df, (2wMG)y df, 

It's important to keep in mind that, because of the dependence of (y 2 ) on fi(y,u;t) (it 
is a functional), this is not simply the equation of a distribution of independent harmonic 
oscillators. 

We have seen that explicit N dependence disappears from the equations of motion in the 
fluid limit. In the continuum limit it's convenient, and in fact necessary, to introduce scaled 
quantities r, h, a, I 2 , 

R h H A 2 L 2 

r ■= -7= , h:= — - , a:= -= , / := — — . (17) 



/~N mN ' m 2 N 2 

Then their limiting values are 

r 2 = (y 2 ) ,h=\ (u 2 ) + ar , a = , I 2 = (y 2 ) (u 2 ) - (yu) 2 . (18) 

A remarkable feature of the fluid limit is that the radial equation is still preserved, i.e. 
substituting the scaled variables gives the identical differential equation in terms of the 
scaled variables, so that r(t) is still a periodic function of time. With a little work, we find 
from the Vlasov equation for // that both h and I 2 are integral invariants of the fluid motion 
defined by Eq. ()15j) . In common with the OGS, there are an infinite number of stationary 
solutions of Eq. (|TH|) . We can construct a family of such solutions by defining the combined 
energy-like function (and functional) g : 

9{y^h\) -.= \u 2 + \^=. (19) 

1 1 v(y) 
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Note that || is the local velocity and is the local acceleration. Therefore any integrable 
functional of g satisfies the Vlasov equation and yields a stationary solution fi[g]. Note also 
that (g) 7^ (h). It is also possible to use the invariance of h and I 2 to construct the family 
of maximum entropy solutions. 

The standard derivation of the Vlasov equation starts with the BBGKY hierarchy for a 
system of pairwise interacting particles and proves that, in the Vlasov limit, the higher order 
distribution functions factor into products of the singlet distribution / jll|, so that pairs 
of particles become statistically independent. In fact, we implicitly made this assumption 
in determining the expression for I 2 . However, in contrast with the OGS, the EIS does not 
have pairwise additive interactions so it may be possible to construct additional continuum 
representations which include correlations. This possibility is suggested by the form of the 
generalized angular momentum, Eq.(j3J), which depends on the coordinates and velocities of 
pairs of particles. Note that, in the continuum limit, both l p := (yu' — y'u) and its square 

(2) 

are invariants of the motion. Let /) (x, u; x', u ; t) be a normalized distribution in the four 
dimensional pair space and assume that it satisfies the extended Vlasov equation 

dfj 2) u djf u , dfP _ (2irMG)y gjf _ {2^MG)y' dff = Q 
Ot U dy U dy' ^W) du y^/ 2 ) du' 
Then, as in the standard formulation, one class of solutions can always be constructed by 
the product fi(y,u;t)fi(y',u';t). However, there are now other possibilities. For example, 
the integrable functionals of l p (e.g. ff^ 1 ~ exp(— cl 2 )) will generate a family of stationary 
solutions in which pairs of points are obviously correlated. It will be interesting to look for 
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connections with earlier studies of correlation in the OGS 

Of course, due to their dynamical instability, nearly all orbits generated by the OGS will 
diverge from their EIS counterparts sooner or later. To get an idea of how well the orbits 
generated by the EIS follow the those of the OGS, we carried out numerical simulations of the 
OGS and integrated the radial equation numerically for the EIS. For the initial conditions we 
selected random initial positions and velocities that are uniformly distributed in a rectangle 
in (x, v). We chose units of mass and time where M := mN = 1 and 2irG = l.We made 
small translations to guarantee that the total momentum and center of mass are null, and 
scaled the positions and velocities so that the total energy is exactly 0.75 and the virial ratio 
(twice the kinetic energy divided by the potential energy) is exactly 1.0. In each case we 
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FIG. 2: Comparison of the radial coordinate R(t) computed for the OGS and EIS with the same 
initial condition. Note the close agreement for this initial condition. 

computed the initial values of R 2 and L 2 which are all that is required to solve the radial 
equation, Eq.ffTT]). 

Below we present graphs comparing the theory, from the EIS, and simulations. Although 
the results are not exhaustive, there seems to be a well established trend. In cases where the 
simulations exhibit a very regular variation of the radius, there is a remarkable agreement 
between theory (as given by the EIS) and simulation. This can be seen in Figure El For 
other, equally likely, initial conditions the simulated R(t) is not regular and, therefore, is 
not well represented by the theory (see Figure EJ). However, even in these cases, we find that 
the maxima and minima of the simulations line up fairly closely with those of the theory, 
i.e. they are in phase. This seems to be a general trend whether we are modelling the initial 
behavior or letting the system anneal for a while before comparisons are made. 

Starting from the N particle one- dimensional self-gravitating system (OGS), we have 
shown how to construct an exactly integrable Hamiltonian system of arbitrary dimension 
(EIS) in which the motion is always confined to a plane in the configuration space. The 
system consists of a single point which experiences an acceleration of constant magnitude 
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FIG. 3: Evolution of the radial coordinate for the OGS from three different initial conditions 
selected using the identical algorithm. 

that is always directed towards a fixed point (say the origin) of iV dimensional Euclidean 
space. The exact solution of the equations of motion can be expressed as a linear combina- 
tion of elliptic integrals. In general, the motion fills an annulus in the plane between two 
turning points, but degenerate cases are possible where the motion is exactly circular (max- 
imum possible angular momentum), resulting in stationary motion, or linear (null angular 
momentum), resulting in the periodic collapse and expansion of the system. By imposing 
the Vlasov limit, we have shown how to construct two types of continuum models for the 
EIS, one of which is capable of representing correlations between particles. Simulations of 
the OGS have been compared with solutions of the EIS with identical initial conditions. For 
cases where the dynamical simulations are regular, the agreement between the two models 
is outstanding. For other less regular behavior, extrema in the simulations seem to occur at 
similar times as the periodic extrema of the integrable model, so the solutions are in phase. 
From the preliminary results obtained so far, it appears that the existence of substructure 
in the gravitational system is responsible for the lack of regularity. Further work concern- 
ing the influence of population, the graininess of distributions, and the theory and possible 
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applications of the continuum representations of the EIS is ongoing. In addition we are 
investigating the extension of the theory to higher dimensional systems. 
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